Complex visualizations

data %>% glimpse

## Rows: 671
## Columns: 26
## $ birth    <dbl> 81.511, 81.514, 81.552, 81.558, 81.593, 81.602, 81.610, NA, N…
## $ exit     <dbl> 81.604, 81.539, 81.552, 81.667, 81.599, 81.771, 81.697, NA, N…
## $ hospstay <int> 34, 9, -2, 40, 2, 62, 32, NA, NA, 28, 38, NA, 62, 69, 1, 93, …
## $ lowph    <dbl> NA, 7.250000, 7.059998, 7.250000, 6.969997, 7.189999, 7.32000…
## $ pltct    <int> 100, 244, 114, 182, 54, NA, 282, NA, NA, 153, 229, NA, 182, 3…
## $ race     <fct> white, white, black, black, black, white, black, NA, NA, blac…
## $ bwt      <int> 1250, 1370, 620, 1480, 925, 940, 1255, 600, 700, 1350, 1310, …
## $ gest     <dbl> 35.0, 32.0, 23.0, 32.0, 28.0, 28.0, 29.5, 26.0, 24.0, 34.0, 3…
## $ inout    <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn      <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 1, 1, 0…
## $ lol      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ magsulf  <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ meth     <int> 0, 1, 0, 1, 0, 1, 1, NA, NA, 1, 0, NA, 0, 0, 1, 1, 1, 1, 1, 1…
## $ toc      <int> 0, 0, 1, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 1, 1, 0, 0, 0, 1…
## $ delivery <fct> abdominal, abdominal, vaginal, vaginal, abdominal, abdominal,…
## $ apg1     <int> 8, 7, 1, 8, 5, 8, 9, NA, NA, 4, 6, NA, 6, 6, 2, 4, 8, 7, 1, 8…
## $ vent     <int> 0, 0, 1, 0, 1, 1, 0, NA, NA, 0, 1, NA, 0, 0, 1, 1, 0, 0, 1, 1…
## $ pneumo   <int> 0, 0, 0, 0, 1, 0, 0, NA, NA, 0, 0, NA, 1, 0, 1, 0, 0, 0, 1, 1…
## $ pda      <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld      <int> 0, 0, NA, 0, 0, 0, 0, NA, NA, 0, 0, NA, 1, 0, 0, 1, 0, 0, 0, …
## $ pvh      <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ivh      <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ipe      <fct> NA, NA, NA, NA, NA, absent, NA, NA, absent, NA, NA, NA, absen…
## $ year     <dbl> 81.51196, 81.51471, 81.55304, 81.55847, 81.59406, 81.60229, 8…
## $ sex      <fct> female, female, female, male, female, female, female, NA, NA,…
## $ dead     <int> 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0…

Задание 1

Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.

cpdata <- data %>%
  select(where(~sum(is.na(.))<=100)) %>% 
  drop_na()

cpdata %>% head

##    birth   exit hospstay    lowph pltct  race  bwt gest        inout twn
## 1 81.514 81.539        9 7.250000   244 white 1370 32.0 born at Duke   0
## 2 81.558 81.667       40 7.250000   182 black 1480 32.0 born at Duke   0
## 3 81.593 81.599        2 6.969997    54 black  925 28.0 born at Duke   0
## 4 81.610 81.697       32 7.320000   282 black 1255 29.5 born at Duke   0
## 5 81.624 81.700       28 7.160000   153 black 1350 34.0 born at Duke   0
## 6 81.626 81.730       38 7.039997   229 white 1310 32.0 born at Duke   0
##    delivery apg1 vent pneumo pda cld     year    sex dead
## 1 abdominal    7    0      0   0   0 81.51471 female    0
## 2   vaginal    8    0      0   0   0 81.55847   male    0
## 3 abdominal    5    1      1   0   0 81.59406 female    1
## 4   vaginal    9    0      0   0   0 81.61053 female    0
## 5 abdominal    4    0      0   0   0 81.62421 female    0
## 6   vaginal    6    1      0   0   0 81.62695   male    0

Задание 2

Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.

В оставшихся данных мы можем увидеть несколько категориальных переменных: twn (multiple gestation), vent (assisted ventilation used), pneumo (pneumothorax occurred), pda (patent ductus arteriosus detected), cld (on suppl. oxygen at 30 days), dead.

str(cpdata)

## 'data.frame':    531 obs. of  19 variables:
##  $ birth   : num  81.5 81.6 81.6 81.6 81.6 ...
##  $ exit    : num  81.5 81.7 81.6 81.7 81.7 ...
##  $ hospstay: int  9 40 2 32 28 38 62 69 1 93 ...
##  $ lowph   : num  7.25 7.25 6.97 7.32 7.16 ...
##  $ pltct   : int  244 182 54 282 153 229 182 361 378 255 ...
##  $ race    : Factor w/ 4 levels "black","native American",..: 4 1 1 1 1 4 1 4 4 1 ...
##  $ bwt     : int  1370 1480 925 1255 1350 1310 1110 1180 970 770 ...
##  $ gest    : num  32 32 28 29.5 34 32 28 28 28 26 ...
##  $ inout   : Factor w/ 2 levels "born at Duke",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ twn     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ delivery: Factor w/ 2 levels "abdominal","vaginal": 1 2 1 2 1 2 2 1 2 2 ...
##  $ apg1    : int  7 8 5 9 4 6 6 6 2 4 ...
##  $ vent    : int  0 0 1 0 0 1 0 0 1 1 ...
##  $ pneumo  : int  0 0 1 0 0 0 1 0 1 0 ...
##  $ pda     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ cld     : int  0 0 0 0 0 0 1 0 0 1 ...
##  $ year    : num  81.5 81.6 81.6 81.6 81.6 ...
##  $ sex     : Factor w/ 2 levels "female","male": 1 2 1 1 1 2 2 2 1 2 ...
##  $ dead    : int  0 0 1 0 0 0 0 0 1 0 ...

cpdata <- cpdata %>%
  mutate_at(vars(twn, vent, pneumo, pda, cld, dead), as.factor)

Строим графики распределений:

distributions <- list()
nums <- cpdata %>% select(is.double | is.integer) %>% colnames()
idx <- 0
for (variable in nums) {
  idx <- idx + 1
  p <- ggplot(cpdata, aes_string(x=variable)) +
    theme_minimal()
  
  p <- if (idx < 3) {
    p + geom_density(aes_string(fill='inout'), alpha=0.5) # ifelse почему то возвращает dataframe, а не ggplot объект. При этом при указании alpha в функции ggplot, параметр не распространяется на добавочные функции geom_density в рамках if else, странно
  } else {
    p + geom_density(fill='darkred', alpha=0.5)
    }
  
  distributions[[idx]] <- p
}

wrap_plots(distributions, axes='collect_y')

Попробуем удалить выбросы и посмотреть как изменятся распределения

remove_outliers <- function(data, column) {
  
  Q1 <- quantile(data[[column]], 0.25)
  Q3 <- quantile(data[[column]], 0.75)
  
  IQR <- Q3 - Q1
  
  lower_bound <- Q1 - 1.5 * IQR
  upper_bound <- Q3 + 1.5 * IQR
  
  filtered_data <- data[data[[column]] >= lower_bound & data[[column]] <= upper_bound, ]
  
  return(filtered_data)
}

cpdata_noout <- cpdata
for (variable in nums) {
  cpdata_noout <- remove_outliers(cpdata_noout, variable)
}

distributions_noout <- list()
nums <- cpdata_noout %>% select(is.double | is.integer) %>% colnames()
idx <- 0
for (variable in nums) {
  idx <- idx + 1
  p <- ggplot(cpdata_noout, aes_string(x=variable)) +
    theme_minimal()
  
  p <- if (idx < 3) {
    p + geom_density(aes_string(fill='inout'), alpha=0.5) # ifelse почему то возвращает dataframe, а не ggplot объект. При этом при указании alpha в функции ggplot, параметр не распространяется на добавочные функции geom_density в рамках if else, странно
  } else {
    p + geom_density(fill='darkred', alpha=0.5)
    }
  
  distributions_noout[[idx]] <- p
}

wrap_plots(distributions_noout, axes='collect_y')

В основном выбросы влияли на переменную hospstay, теперь она выглядит намного лучше(отсутствуют отрицательные значений(скорее всего ошибка ввода), график плотности стал менее вытянутым). Будем работать с данными без выбросов.

Задание 3

Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?

set.seed(123)

shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == 'born at Duke'])

## # A tibble: 1 × 3
##   variable                                                     statistic p.value
##   <chr>                                                            <dbl>   <dbl>
## 1 "cpdata_noout$lowph[cpdata_noout$inout == \"born at Duke\"]"     0.984 1.45e-4

shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == "transported"])

## # A tibble: 1 × 3
##   variable                                                    statistic p.value
##   <chr>                                                           <dbl>   <dbl>
## 1 "cpdata_noout$lowph[cpdata_noout$inout == \"transported\"]"     0.967  0.0562

if (shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == 'born at Duke'])$p.value > 0.05 && 
    shapiro_test(cpdata_noout$lowph[cpdata_noout$inout == "transported"])$p.value > 0.05) {
  # Если данные нормально распределены, используем t-тест
  test_result <- t_test(lowph ~ inout, data = cpdata_noout)
  print('Использован Т-тест')
} else {
  # Если данные не нормально распределены, используем U-тест Манна-Уитни
  test_result <- wilcox_test(lowph ~ inout, data = cpdata_noout)
  print('Использован тест Манна-Уитни')
}

## [1] "Использован тест Манна-Уитни"

# Выводим результаты теста
print(test_result)

## # A tibble: 1 × 7
##   .y.   group1       group2         n1    n2 statistic           p
## * <chr> <chr>        <chr>       <int> <int>     <dbl>       <dbl>
## 1 lowph born at Duke transported   413    73     20657 0.000000446

ggplot(cpdata_noout, aes(x = inout, y = lowph)) +
  geom_boxplot() +
  stat_compare_means(method = "wilcox.test") +
  labs(title = "Сравнение значений lowph между группами",
       x = "Группы",
       y = "Значения lowph")

Интерпретация: тк среднее значение low pH статистически значимо ниже у транспортированных пациентов и его снижение положительно ассоциированно с риском смерти, то можно предположить, что транспортировка (условия транспортировки и/или время затраченное на нее) ухудшают состояние пациента и снижают его шансы на выживание.

Задание 4

Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

ftask <- cpdata_noout %>% select(-c(birth, year, exit)) %>% select(is.double|is.factor)

cormat <- ftask %>% mutate(across(where(is.factor), as.numeric)) %>% cor()

corrplot(cormat)

ggcorrplot(cormat)

Задание 5

Постройте иерархическую кластеризацию на этом датафрейме.

ftask <- cpdata_noout %>% select(-c(birth, year, exit)) %>% select(is.double|is.factor) %>% mutate(across(where(is.factor), as.numeric))

distmat <- dist(t(scale(ftask)), method = "euclidean")

hc <- hclust(distmat, method = "complete")

plot(hc, main = "Hierarchical Clustering Dendrogram", xlab = "Variables", sub = "", cex = 0.8)

pheatmap(
  scale(ftask),
  clustering_distance_rows = "euclidean",
  clustering_distance_cols = "euclidean",
  clustering_method = "complete",
  main = "Heatmap with Hierarchical Clustering"
)

Интерпретация:

Группа умерших людей подсвечена красным в колонке dead. По хитмапу можно сказать, что pda, pneumo, vent имеют хорошие перекрытия с dead (и являются частью одной клады), при этом lowph тоже имеют хорошую отрицательную корреляцию с исходом, эти переменные можно было бы попробовать использовать в качестве предикторов.

Задание 7

Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?

Переменные в нашем датасете имеют разные диапозоны значений, так что шкалирование необходимо, чтобы исключить доминирование “больших” переменных над малыми

df_scaled <- scale(ftask)
pca_result <- prcomp(df_scaled, center = TRUE, scale. = TRUE)
summary(pca_result)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.7509 1.1969 1.05804 1.04516 0.97608 0.93719 0.90775
## Proportion of Variance 0.2555 0.1194 0.09329 0.09103 0.07939 0.07319 0.06867
## Cumulative Proportion  0.2555 0.3749 0.46815 0.55918 0.63857 0.71176 0.78043
##                            PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.83482 0.76392 0.74432 0.66321 0.60038
## Proportion of Variance 0.05808 0.04863 0.04617 0.03665 0.03004
## Cumulative Proportion  0.83851 0.88714 0.93331 0.96996 1.00000

Видно, что для объяснения 75 процентов дисперсии необходимо взять целых 7(из 12) компонент, значит, что в данных нету переменных, которые бы однозначно их описывали.

Задание 8

# Выполнение PCA
pca_result <- prcomp(df_scaled, center = TRUE, scale. = TRUE)
# Построение biplot
p <- fviz_pca_biplot(
  pca_result,
  text = paste0('Patient ID: ', rownames(ftask)),
  geom.ind = "point",
  palette = c("blue", "red"),
  repel = TRUE
) + labs(title = "PCA Biplot") +
  geom_point(aes(col= as.factor(ftask$dead)), text = paste0('Patient ID: ', rownames(ftask))) +
  scale_color_discrete('Patient', labels = c('1' = 'Alive', '2' = 'Dead'))

p

Задание 9

Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.

plot_ly(
    x = pca_result$x[,1],
    y = pca_result$x[,2],
    text = pca_result$x %>% rownames(),
    hovertemplate = paste('<br><b>Patient ID</b>: <b>%{text}</b>'),
    hoverinfo = 'text',
    type = 'scatter',
    mode = 'markers',
    color = fct_recode(as.factor(ftask$dead), 'Alive' = '1', 'Dead' = '2')
  ) %>% 
  layout(title = "PCA Biplot",
         xaxis = list(title = "Dim2(11%)"),
         yaxis = list(title = "Dim1(23.6%)"),
         legend = list(title = list(text = "Patient Status")))